Setting the data
# Get palestinian cities, convert it to sf dataframe, assign projection to the coordinates and save it
data(world.cities)
palestine_cities <- world.cities %>% filter(country.etc == "Palestine")
class(palestine_cities)
## [1] "data.frame"
palestine_cities <- as_tibble(palestine_cities)
write_excel_csv(palestine_cities, file = "D:/R_Gaza/pal_map/data/palestinians_cities.csv")
names(palestine_cities) #get the names for the coordinate variables
## [1] "name" "country.etc" "pop" "lat" "long"
## [6] "capital"
palestine_cities_sf <- st_as_sf(palestine_cities, coords = c("lat", "long")) #convert to sf df tibble
class(palestine_cities_sf) #check if class has changed to sf dataframe tibble
## [1] "sf" "tbl_df" "tbl" "data.frame"
st_crs(palestine_cities_sf) #check for geographic projections
## Coordinate Reference System: NA
st_crs(palestine_cities_sf) <- 4326 #assign geographic projection (4326, i.e., longlat WGS84)
st_crs(palestine_cities_sf) #check if the geographic projection was assigned
## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
st_crs(palestine_cities_sf)$proj4string
## [1] "+proj=longlat +datum=WGS84 +no_defs"
palestine_cities_sf #check if "lat" and "long" variables were changed to "geometry" variable.
## Simple feature collection with 244 features and 4 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 31.27 ymin: 34.25 xmax: 32.55 ymax: 35.45
## Geodetic CRS: WGS 84
## # A tibble: 244 × 5
## name country.etc pop capital geometry
## * <chr> <chr> <int> <int> <POINT [°]>
## 1 'Abasan al-Jadidah Palestine 5629 0 (31.31 34.34)
## 2 'Abasan al-Kabirah Palestine 18999 0 (31.32 34.35)
## 3 'Abud Palestine 2456 0 (32.03 35.07)
## 4 'Abwein Palestine 3434 0 (32.03 35.2)
## 5 'Ajjah Palestine 5143 0 (32.37 35.2)
## 6 'Anabta Palestine 7311 0 (32.31 35.12)
## 7 'Anata Palestine 9149 0 (31.81 35.28)
## 8 'Anin Palestine 3717 0 (32.5 35.17)
## 9 'Anzah Palestine 2006 0 (32.37 35.22)
## 10 'Aqbah Jabbar Palestine 6337 0 (31.83 35.45)
## # ℹ 234 more rows
st_write(palestine_cities_sf, "D:/R_Gaza/pal_map/data/palestinian_cities",
driver = "ESRI Shapefile", delete_layer = TRUE) # save the sf data frame as shapefile
## Deleting layer `palestinian_cities' using driver `ESRI Shapefile'
## Writing layer `palestinian_cities' to data source
## `D:/R_Gaza/pal_map/data/palestinian_cities' using driver `ESRI Shapefile'
## Writing 244 features with 4 fields and geometry type Point.
unzip(zipfile = "D:/R_Gaza/gaza/data/localities_palestine_sf.zip" , exdir ="data/localities_palestine")
localities_palestine_sf <- st_read("data/localities_palestine")
## Multiple layers are present in data source D:\R_Gaza\pal_map\data\localities_palestine, reading layer `%D8%AA%D8%AC%D9%85%D8%B9%D8%A7%D8%AA_%D9%81%D9%84%D8%B3%D8%B7%D9%8A%D9%86___Localities_Palestine'.
## Use `st_layers' to list all layer names and their type in a data source.
## Set the `layer' argument in `st_read' to read a particular layer.
## Reading layer `%D8%AA%D8%AC%D9%85%D8%B9%D8%A7%D8%AA_%D9%81%D9%84%D8%B3%D8%B7%D9%8A%D9%86___Localities_Palestine' from data source `D:\R_Gaza\pal_map\data\localities_palestine' using driver `ESRI Shapefile'
## Simple feature collection with 619 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 3809204 ymin: 3661924 xmax: 3960023 ymax: 3836009
## Projected CRS: WGS 84 / Pseudo-Mercator
localities_palestine_sf
## Simple feature collection with 619 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 3809204 ymin: 3661924 xmax: 3960023 ymax: 3836009
## Projected CRS: WGS 84 / Pseudo-Mercator
## First 10 features:
## OBJECTID REGIONCODE DISTCODE GOVCODE LOCCODE NAMEAR
## 1 1 2 24 2470 24703430 عبسان الجديدة
## 2 2 2 24 2470 24703470 خزاعة
## 3 3 2 24 2465 24653200 مخيم دير البلح
## 4 4 2 24 2465 24653275 وادي السلقا
## 5 5 2 24 2455 24552681 أم النصر (القرية البدوية)
## 6 6 2 24 2465 24653065 مخيم النصيرات
## 7 7 2 24 2470 24703420 خانيونس
## 8 8 2 24 2470 24703425 بني سهيلا
## 9 9 2 24 2470 24703445 عبسان الكبيرة
## 10 10 2 24 2470 24703485 الفخاري
## NAMEEN
## 1 'Abasan al Jadida
## 2 Khuza'a
## 3 Deir al Balah Camp
## 4 Wadi as Salqa
## 5 Um Al-Nnaser (Al Qaraya al Badawiya)
## 6 An Nuseirat Camp
## 7 Khan Yunis
## 8 Bani Suheila
## 9 'Abasan al Kabira
## 10 Al Fukhkhari
## NOTES
## 1 <NA>
## 2 <NA>
## 3 <NA>
## 4 <NA>
## 5 <NA>
## 6 <NA>
## 7 تشمل المَواصِي (خانيونس)، وقِيزَان النَجَّار، وقَاعْ الخَرَابَة، وقَاعْ القُرِين، وأُمُّ كَمِيل، وأُّمُّ الكِلاب
## 8 <NA>
## 9 <NA>
## 10 <NA>
## Shape__Are Shape__Len geometry
## 1 4919436.7 11010.191 MULTIPOLYGON (((3824072 367...
## 2 9193057.5 14360.481 MULTIPOLYGON (((3826180 367...
## 3 337604.1 3122.714 MULTIPOLYGON (((3823061 368...
## 4 8752155.1 12611.773 MULTIPOLYGON (((3827289 368...
## 5 4306257.9 9547.706 MULTIPOLYGON (((3844490 370...
## 6 1099596.7 10312.221 MULTIPOLYGON (((3827522 369...
## 7 72893548.1 53756.789 MULTIPOLYGON (((3821981 367...
## 8 9181044.3 16897.216 MULTIPOLYGON (((3823716 367...
## 9 17770337.9 20809.409 MULTIPOLYGON (((3823887 367...
## 10 12900017.9 17496.895 MULTIPOLYGON (((3820790 367...
class(localities_palestine_sf)
## [1] "sf" "data.frame"
st_crs(localities_palestine_sf)$proj4string
## [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"
localities_palestine_WGS84 <- st_transform(localities_palestine_sf, crs = 4326)
st_crs(localities_palestine_WGS84)$proj4string
## [1] "+proj=longlat +datum=WGS84 +no_defs"
st_write(localities_palestine_WGS84, "D:/R_Gaza/pal_map/data/localities_palestine",
driver = "ESRI Shapefile", delete_layer = TRUE)
## Deleting layer `localities_palestine' using driver `ESRI Shapefile'
## Writing layer `localities_palestine' to data source
## `D:/R_Gaza/pal_map/data/localities_palestine' using driver `ESRI Shapefile'
## Writing 619 features with 10 fields and geometry type Multi Polygon.
st_crs(localities_palestine_WGS84)$proj4string
## [1] "+proj=longlat +datum=WGS84 +no_defs"
st_crs(palestine_cities_sf)$proj4string
## [1] "+proj=longlat +datum=WGS84 +no_defs"
Creating maps
ggplot(localities_palestine_WGS84)+
geom_sf(colour = "black")

col_fun <- colorFactor(c("#337bd7", "#c8143b"), domain = NULL)
itamar <- st_point(c(34.82375, 31.89530))
itamarpopup <- c("Where I.J.R was born")
pal_popup <- paste0("<strong>City Name:</strong>", palestine_cities$name)
palestine_map <- leaflet(localities_palestine_WGS84) %>%
addPolygons(stroke = TRUE, weight = 0.7, color = "black",
fillColor = ~col_fun(REGIONCODE),
group = "Regions", fillOpacity = 0.8) %>%
addCircleMarkers(lng = 34.82375 , lat = 31.89530, group = "My House",
popup = itamarpopup) %>%
addCircleMarkers(data = palestine_cities, group = "Cities", popup = pal_popup,
radius = 0.01, color = "#2c120c") %>%
addTiles(group = "OSM") %>%
addProviderTiles(provider = "Stadia.OSMBright", group = "OSMBright") %>%
addProviderTiles(provider = "Stadia", group = "Stadia") %>%
addProviderTiles(provider = "Stadia.Outdoors", group = "Outdoors") %>%
addProviderTiles(provider = "Stadia.StamenTonerLite", group = "TonerLite") %>%
addProviderTiles(provider = "CartoDB", group = "CartoDB") %>%
addProviderTiles(provider = "CartoDB.Voyager", group = "CartoDB.Voyager") %>%
addProviderTiles(provider = "Esri", group = "Esri") %>%
addLayersControl(baseGroups = c("OSM", "OSMBright", "Stadia", "Outdoors", "TonerLite",
"CartoDB", "CartoDB.Voyager", "Esri"),
overlayGroups = c("Regions", "Where I Live", "Cities"))
## Assuming "long" and "lat" are longitude and latitude, respectively
palestine_map